Force-matched embedded-atom method potential for niobium 
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Large-scale simulations of plastic deformation and phase transformations in alloys require reliable 
classical interatomic potentials. We construct an embedded-atom method potential for niobium as 
the first step in alloy potential development. Optimization of the potential parameters to a well- 
converged set of density-functional theory (DFT) forces, energies, and stresses produces a reliable 
and transferable potential for molecular dynamics simulations. The potential accurately describes 
properties related to the fitting data, and also produces excellent results for quantities outside the 
fitting range. Structural and elastic properties, defect energetics, and thermal behavior compare 
well with DFT results and experimental data, e.g., DFT surface energies are reproduced with less 
than 4% error, generalized stacking-fault energies differ from DFT values by less than 15%, and the 
melting temperature is within 2% of the experimental value. 

PACS numbers: 34.20.Cf, 62.20.-x, 65.40.-b, 61.72.J- 
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I. INTRODUCTION 

Niobium's low density, high melting temperature, and 
biocompatibility make it an attractive material for al- 
loy design. Nb alloys are promising candidate materials 
for a wide variety technological applications. Multifunc- 
tional Ti-based "gum metal" alloys with substantial Nb 
concentrations exhibit remarkable properties and unique 
deformation behavior—. Attempts to increase the oper- 
ating temperatures, and hence efficiencies, of turbine en- 
gines have prompted interest in designing Nb-based su- 
peralloys 2 -. Non-toxic Ti-Nb based shape-memory alloys 
offer an alternative to Ti-Ni alloys for biomedical applica- 
tions 3 . Accurate and efficient computational models will 
aide in the microscopic understanding of deformation and 
transformation processes in these materials. 

Advances in computational hardware and algorithms 
allow application of first-principles methods to increas- 
ingly complex problems. However, there remain calcula- 
tions beyond the realm of ah initio methods. Meaning- 
ful simulations of processes involving long-ranged strain 
fields or long-wavelength fluctuations require large num- 
bers of atoms. Therefore, methods must be developed 
that scale more favorably with system size than first- 
principles calculations, while retaining a high degree of 
accuracy. The computational cost for simulations based 
on short-ranged classical potentials scales linearly with 
system size, allowing routine simulations of millions of 
atoms. However, the potentials must be carefully con- 
structed and thoroughly tested to ensure that they yield 
reliable results. We construct a classical potential for 
large-scale Nb simulations, and subsequent incorporation 
into potentials for alloys including Ti-Nb. 

A number of authors have developed classical Nb po- 
tentials based on analytic functions^—. Analytic poten- 
tials are typically fit to experimental data for a small 
number of properties. These potentials reproduce the 
fit data with high accuracy, but they often have lim- 
ited transferability and can produce inaccurate forces 



for molecular dynamics (MD) simulations. The force- 
matching method proposed by Ercolessi and Adamsi^ of- 
fers an alternative way to construct potentials. The func- 
tions are parameterized by cubic splines, and the spline 
knots are fit to a large number of forces from density- 
functional theory (DFT) 17 i 18 calculations and usually ex- 
perimental data as well. Including force data from dif- 
ferent thermodynamic conditions improves accuracy and 
transferability for a larger range of simulations. 

We use the force-matching method to develop a cu- 
bic spline-based embedded-atom method (EAM) poten- 
tia l 19 ' 20 for Nb. The potential is fit to a database of 
DFT forces, energies, and stresses from ah initio molec- 
ular dynamics (MD) simulations. We do not fit to any 
experimental data, since it may be inconsistent with the 
DFT information. Instead, we use experimental data and 
DFT results to test the accuracy of the potential. Section 
II discusses our DFT database calculations and poten- 
tial optimization process. We utilize the force-matching 
program potfi t 21 ' 22 to optimize the spline knots to the 
DFT database. Section III presents EAM calculations of 
structural and elastic properties, defects, and thermal be- 
havior. We compare the results to DFT calculations and 
experimental data. These calculations demonstrate the 
potential's ability to describe properties related to the 
fitting data, as well as transferability to behavior beyond 
the fitting range. 



II. OPTIMIZATION OF THE 
EMBEDDED-ATOM METHOD POTENTIAL 

A. Embedded-atom method interatomic potentials 

EAM potential s 19 ' 20 overcome limitations associated 
with simple pairwise interatomic potentials in simula- 
tions of metallic systems. Pair potentials yield a num- 
ber of incorrect predictions for transition metals^ 3 -, in- 
cluding bond energies that are independent of the lo- 
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cal bonding environment, a zero value for the Cauchy 
pressure (Ci 2 — C44 = 0), and the equivalence of the 
cohesive energy with the unrelaxed vacancy formation 
energy. Supplementing the pairwise interaction with a 
volume-dependent term removes some of these undesir- 
able feature a 23 ' 24 , but the volume is ill-defined near de- 
fects such as cracks and surfaces. EAM potentials over- 
come these difficulties by implicitly including many-body 
interactions, and requiring the local "electronic density" 
as input rather than volume. 

The EAM formalism is based on ideas from DF T 17 ' 18 . 
The energy required to embed an impurity atom Z in 
a solid at position R is a unique functional -Fz,R,[n] of 
the electronic density n of the solid before the impurity 
is adde d 25 ' 26 . The embedded-atom method views each 
atom in the solid as an impurity embedded in a host solid 
made up of the remaining atoms. The energy functional 
is approximated by a potential energy function with two 
terms: (1) a sum of pairwise interactions <j) SiSj (\ri — rj\) 
between atoms i and j, and (2) a sum of embedding en- 
ergies F Si (rij) for each atom i that depend on the local 
electronic density rii the atom sees from its surrounding 
neighbors. This local electronic density is a sum of radi- 
ally symmetric electronic-density functions p Sj (\ri — r j- 1 ) 
arising from the atoms j surrounding a given atom z, 

ni = ^2p Si {rij), (1) 

where = |r; — rj \ is the distance between atoms i and 
j. The total potential energy is 

£ = X)W r «) ( 2 ) 

i<j i 

where the subscripts &j and Sj indicate that the functions 
depend on the species of the atoms. 

Equations |T]) and are general and hold for multi- 
component systems. The energy expression simplifies for 
monatomic systems, 

i<j i 

where 

n i = ^2,p{r ij ). (4) 

Thus, for a single component system the three func- 
tions (j), F, and p must be determined (whereas two- 
component alloys require seven functions). EAM poten- 
tials are implemented in many freely-available MD codes, 
e.g., imd2L2S, lammps22, and OHMmA We perform our 
EAM calculations using all three of these codes. 



B. Database of DFT forces, energies, and stresses 



We use the force-matching method of Ercolessi and 
Adam&i 6 - to obtain accurate potentials for molecular dy- 
namics simulations. Force-matched potentials are fit to 
forces from DFT calculations and typically physical prop- 
erties from either DFT calculations or experiment. Here 
we include only DFT data in our fitting database to avoid 
conflicting information. We use the force-matching pro- 
gram potfi t 21 i 22 to optimize the EAM functions to a 
database of DFT forces, energies per atom, and stresses 
for Nb from the configurations listed in Table HI Fitting 
to DFT data from configurations under different temper- 
ature and strain conditions extends the applicability of 
the potential to a wide range of simulations. 

The DFT calculations are performed using the plane- 
wave program vasp- 3 ^. The Perdew-Burke-Ernzerhof 
(PBE) generalized-gradient approximation (GGA) func- 
tional' 52 accounts for the electronic exchange-correlation 
energy, and a projector augmented-wave (PAW) pseu- 
dopotential— generated by Kresse^i represents the nu- 
cleus and core electrons. Along with the five valence 
states, the 4s and 4p semicore states are treated explicitly 
to accurately describe interactions at small interatomic 
separations. The elastic constants also agree better with 
experiment when more electronic states are included. 

The database is calculated in two steps. First, ab ini- 
tio (MD) simulations generate realistic atomic trajecto- 
ries for various temperature and strain conditions. The 
simulation supercells contain 124-128 atoms, depending 
on the structure. These calculations use a relatively low 
convergence criteria to reduce the computational burden. 
A single k point is used, and the plane-wave cutoff en- 
ergy is set to the default value of 219.927 eV from the 
VASP pseudopotential file. Order-one Methfessel-Paxton 
smearing^ 5 , is used with a smearing width of 0.10 eV. The 
MD simulations run for 400 steps with a 3 fs time-step. 

Second, well-converged calculations determine the 
forces, energies per atom, and stresses for the atomic con- 
figurations resulting from the final step of the MD simula- 
tions. We use T-centered fc-point meshes with 30 x 30 x 30 
points per primitive cell and increase the plane- wave cut- 
off energy to 550 eV. The value of the smearing parameter 
is unchanged. The energies are converged to less than 1 
meV/atom. The fitting database contains 1,895 forces 
(5,685 force components), 15 energies per atom, and 90 
stress tensor components from these calculations. Ta- 
ble U lists the configurations in the database, along with 
the weighted relative rms deviations of the EAM force 
magnitudes from the DFT values and the weighted aver- 
age angular deviations of the EAM force directions from 
the DFT force directions (these errors are discussed in 
Sec. EH . 
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TABLE I: Configurations in the fitting database. The "structure" column lists the crystal structure of each configuration. 
The "primitive" and "conventional" labels indicate if the supercell is based on the one-atom primitive cell or the two-atom 
conventional cell. The liquid configuration starts as a bcc lattice and then melts during the ab initio MD simulation. The 
"-^Vatoms" column lists the number of atoms in each configuration. The "T" column lists the temperatures of the ab initio 
MD simulations used to generate the configurations. The "V/Vo" column lists the ratio of the supercell volume to the zero- 
temperature, equilibrium DFT volume of the bcc, fee, or hep structure, respectively. For configuration 13, Vq is the equilibrium 
volume of the bcc structure. The "strain" column indicates the strain applied to the supercells, where M denotes a volume- 
conserving monoclinic strain and O denotes a volume-conserving orthorhombic strain. The "rms deviation" column lists the 
weighted relative rms deviations of the EAM force magnitudes from the DFT values. The "# a v g " column lists the weighted 
average angular deviation of the EAM force directions from the DFT force directions. The weighted relative rms force- magnitude 
deviation and weighted average angular deviation of the forces for all the configurations are 25% and 15.1°, respectively. 



Configuration 


Structure 


A^ atoms 


T(K) 


V/V 


Strain 


rms deviation (%) 


6»avg (deg.) 


1 


bcc, primitive 


125 


300 


0.90 


None 


18 


11.7 


o 
1 


bcc, primitive 




oUU 


i on 


None 


on 


1 O 9 

Iz.o 


3 


bcc, primitive 


125 


300 


1.10 


None 


29 


18.4 


4 


bcc, primitive, vacancy 


124 


300 


1.00 


None 


38 


17.9 


5 


bcc, conventional 


128 


300 


1.00 


2%, M 


20 


14.3 


6 


bcc, conventional 


128 


300 


1.00 


1%, M 


22 


16.1 


7 


bcc, conventional 


128 


300 


1.00 


-1%, M 


21 


14.0 


8 


bcc, conventional 


128 


300 


1.00 


-2%, M 


23 


15.0 


9 


bcc, conventional 


128 


300 


1.00 


2%, O 


23 


16.0 


10 


bcc, conventional 


128 


300 


1.00 


-2%, O 


21 


15.9 


11 


bcc, primitive 


125 


1,200 


1.00 


None 


21 


13.2 


12 


bcc, primitive 


125 


2,200 


1.00 


None 


19 


12.4 


13 


liquid, primitive 


125 


5,000 


1.00 


None 


21 


11.9 


14 


fee, primitive 


125 


300 


1.00 


None 


36 


29.0 


15 


hep, conventional 


128 


300 


1.00 


None 


51 


37.4 



C. Optimization of EAM functions to DFT data 



Generating accurate potentials using the force- 
matching method is an optimization problem in a high- 
dimensional space. The EAM functions are parameter- 
ized by cubic splines, and the program POTFi T 21 i 22 op- 
timizes the spline knots using a combination of simu- 
lated annealing and conjugate-gradientlike algorithms. 
Our potential construction procedure proceeds itera- 
tively. We generate a database of DFT calculations and 
choose an initial set of spline knots. We also specify 
the cutoffs for the functions and the fitting weights for 
the values in the database. Then the optimization algo- 
rithms in potfit adjust the spline knots to minimize the 
weighted error between the database values and the cor- 
responding values produced by the EAM potential. If the 
fitting errors are too large and the potential fails to pro- 
duce satisfactory results for physical properties, we add 
or remove configurations from the database, change the 
fitting weights and cutoffs, and/or modify the number 
and initial values of the spline knots, and refit the poten- 
tial. This optimization and testing process is repeated 
until we obtain accurate potentials. 

In potfit, the fitting error is defined through a least- 
squares target function formed from the differences be- 
tween the EAM and DFT values: 



Z = Zp + Zq, 



(5) 



where 



N A 3 ( pVAM _ pUFT\ ' 



and 



N c ^EAM _ 



A 



DFT\ 



(7) 



Equation (j6|) is the relative deviation of the EAM forces 
from the DFT forces, where 2Va is the number of atoms 
in the fitting database, Fi iX . is the Xjth component of the 
force on atom i, Wi is the weight associated with each 
force, and et is a small number that prevents overweight- 
ing of very small, inaccurate forces. Equation ([7]) is the 
relative deviation between the EAM energies per atom 
and stresses and the DFT values, where Nq is the number 
of energies per atom and stress tensor components, is 
an energy or stress value, Wi is the associated weight, and 
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ti is a small number that prevents overweighting of nu- 
merically small data. The optimal spline knots minimize 
Z. See Brommer and G able r 21 ' 22 for potfit details. 

We determine if there are more parameters in the po- 
tential than the fitting database can support, i.e., over- 
fitting, by calculating the errors for a testing database 
of DFT forces, energies per atom, and stresses for bcc, 
fee, and hep configurations not included in the fitting 
database. If the errors for the testing database are much 
larger than the errors for the fitting database, there are 
likely too many parameters specifying the EAM func- 
tion o 22 ' 36 ' 37 . We also test the optimized potential's abil- 
ity to predict physical properties (see Sec. If the 
potential fails to adequately describe the databases and 
desired properties, we add or remove configurations from 
the databases, modify the fitting weights and cutoffs, 
and/or change the number and initial values of the spline 
knots, and the optimization and testing process repeats. 

Typical of simulated annealing methods, several hun- 
dred iterations of this procedure were required to find a 
small number of reasonable potentials. For the fitting 
database listed in Table U we find that potentials with 
15-20 spline knots for <f> and p and 5-10 spline knots for 
F, and a cutoff radius of 4.75 A for 4> and p produce the 
most physically reasonable results, while giving similar 
force-matching errors of 20-30% for the fitting and test- 
ing databases. The cutoff radius includes first, second, 
and third nearest-neighbor interactions in bcc Nb. The 
cutoffs for F are updated automatically by potfit as p 
changes during optimization. 



D. Optimized potential 

Figure [T] shows the optimized cubic splines of the best 
Nb potential that we found. The pair potential cj> and 
electronic density p are parameterized using seventeen 
equally spaced spline knots while eight spline knots are 
used for the embedding function F. The outer cutoff 
distance for <fi and p is 4.75 A. The shortest interatomic 
distance in the fitting database, which is 2.073 A, de- 
termines the inner cutoff distance. The inner and outer 
cutoffs of F are 0.0775 and 1.000, respectively. 

The function <j> has the expected characteristics for a 
pair potential. The interaction is attractive for large in- 
teratomic separations, and highly repulsive when atoms 
approach too closely. The minimum value for cj> occurs 
at r = 3.197 A. The electronic density p is large for 
small interatomic separations and decreases for r-values 
up to about 3.25 A. Beyond this distance, p ripples and 
decays to zero at 4.75 A. The zero-temperature equilib- 
rium value of n is 0.263. We find that the non-monotonic 
character of p is required for an accurate description of 
Nb. Potentials with smoother p functions found dur- 
ing the optimization and testing procedure yield poor 
results for many properties. The embedding function F 
has positive curvature over most of its range, but there 
is a small region of negative-curvature around n = 0.74. 



This behavior is not ideal but atoms rarely sample ri- 
valries greater than 0.6 even at large temperatures and 
pressures. 

In addition to specifying the spline knots and requir- 
ing continuity of the first and second derivatives of the 
functions at the knots, two boundary conditions must be 
applied to each function to determine all the cubic poly- 
nomial coefficients. The natural boundary condition, i.e., 
a vanishing second derivative, is applied at the inner cut- 
off radius of <fi and p, and at the inner and outer cutoffs 
of F. The remaining boundary conditions are the first 
derivatives of <fi and p are zero at the outer cutoff radius. 
Appendix discusses modifications to <fi, p, and F for small 
and large values of their arguments. These modifications 
improve the performance of the potential at large tem- 
peratures and pressures. Table |TT] lists the spline knots 
and boundary conditions for (f>, p, and F. The potential 
functions are available in tabulated form upon request. 



E. Fitting errors 

The fitting database contains DFT forces, energies, 
and stresses for the configurations listed in Table fl] The 
EAM potential computes the same set of quantities for 
the fixed atomic positions of each configuration and we 
evaluate the deviations of the EAM values from the corre- 
sponding DFT values. The errors associated with numer- 
ically small DFT data are typically much greater than the 
errors from larger data. This is illustrated in Figs. HJa) 
and[2Ib), which show the relative force- magnitude devi- 
ation and angular deviation of each of the EAM forces 
from the DFT database values, versus the DFT force 
magnitudes. Since very small values are inherently inac- 
curate, we weight the terms in the error calculations by 
the magnitudes of the DFT values. The weighted relative 
rms deviation of the energies per atom, stresses, or force 
magnitudes is 



\ 



£< 

i=l 



Q 



EAM 



Q 



DFT 



DFT 



X 100%, (8) 



where Qi is an energy per atom, a stress tensor com- 
ponent, or a force magnitude, and Nq is the respective 
number of such quantities in the database. The scaled 
magnitudes of the DFT data Wi weight the terms in the 
sum, 



UJ, 





Qi 


DFT 




sr N Q 
2^=1 


Qj DFT | 



(9) 



The EAM potential reproduces the energies per atom 
of the configurations with a weighted rms deviation of 
only 0.1%. The diagonal components of the stress ten- 
sors are also accurately reproduced with a 6% weighted 
rms deviation. In contrast, the error for the off-diagonal 
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FIG. 1: The three cubic splines of the EAM potential. The points are optimized spline knots and the solid lines are cubic 
polynomials that interpolate between the knots, (a) The pair potential (p and (b) the "electronic density" p are functions of the 
distance r between pairs of atoms. Both of these functions have seventeen optimized spline knots and a cutoff radius of 4.75 
A, which includes first-, second-, and third-nearest-neighbor interactions in bcc Nb. (c) The embedding function F depends on 
the local electronic density n. Eight optimized spline knots parameterize F. 



TABLE II: The cubic spline knots and boundary conditions. Spline knots 1-17 for <f) and p, and spline knots 1-8 for F are 
optimized by potfit. The adjusted values of knot for <j>, p, and F, and knot 9 for F are also listed (see Appendix). The 
coefficients of the cubic polynomials that interpolate between the knots are determined by requiring continuity of the functions 
and their first and second derivatives, along with the boundary conditions. 



i 


n (A) 


<t>{n) (eV) 


p(n) 


rii 


F(th) (eV) 





1.7383750 


5.644808063640994 


0.683176019233847 


0.000000000000000 


0.000000000000000 


1 


2.0730000 


1.952032491449762 


0.418661384304128 


0.077492938439077 


-3.347285692522362 


2 


2.2403125 


1.094035979464646 


0.248142385672424 


0.209279661519209 


-4.546334492745762 


3 


2.4076250 


0.510885854762808 


0.135151131573890 


0.341066384599341 


-4.893456225397550 


4 


2.5749375 


0.082343335887366 


0.067802030440920 


0.472853107679473 


-4.950236437181159 


5 


2.7422500 


-0.177550651790219 


0.037078599738033 


0.604639830759605 


-4.944970691786193 


6 


2.9095625 


-0.311331931736446 


0.023834158891363 


0.736426553839736 


-4.845699482076931 


7 


3.0768750 


-0.390004380615947 


0.013226669087316 


0.868213276919868 


-4.743717588952880 


8 


3.2441875 


-0.405551151985570 


0.008594239037838 


1.000000000000000 


-4.556142211118433 


9 


3.4115000 


-0.351882201216042 


0.009026077313542 


1.263573446160264 


4.828348385154062 


10 


3.5788125 


-0.251634925091355 


0.013228711231271 






11 


3.7461250 


-0.145378019920633 


0.016102598867695 






12 


3.9134375 


-0.078119761728408 


0.011199412726043 






13 


4.0807500 


-0.047220500113419 


0.007407238328861 






14 


4.2480625 


-0.032830828537903 


-0.002416625008422 






15 


4.4153750 


-0.021236531023427 


-0.002572474995293 






16 


4.5826875 


-0.006495370564318 


-0.000515878027624 






17 


4.7500000 


0.000000000000000 


0.000000000000000 










Boundary conditions 










<t>"(n) = 


4>'(rir) = 










p"(n) = 


p'(r u ) = 










F"(m) = o 


F"(n 8 ) = 









components of the stress tensors is very large. The 
weighted rms deviation for these quantities is 307%. The 
off-diagonal values are very small however, even for the 
strained supercells. Increasing the strain yields larger 
stress values, but the DFT stress-strain curves for Nb 
become non-linear for strains larger than about 2.5%. 



Table U lists the weighted relative rms deviations of the 
force magnitudes for each configuration in the database. 
The weighted relative rms deviation for all the configura- 
tions is 25%. This force-magnitude error is lower than the 
32% error of Li et al.'s force-matched tantalum EAM po- 
tential, and similar to the error of Hennig et al.'s force- 
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F DFT (eV/A) F DFT (eV/A) 6 (degrees) 

FIG. 2: Errors between the 1,895 DFT forces in the fitting database and the corresponding EAM forces. The relative errors in 
the forces tend to decrease as the force magnitudes increase, (a) The percent relative deviation of the EAM force magnitudes 
from the DFT values, versus the magnitudes of the DFT forces. The weighted relative rms deviation of all the forces is 25%. 
The inset shows the absolute errors in the forces, (b) The angular deviations of the EAM force directions from the corresponding 
DFT force directions. The weighted average of all the angular deviations is 15.1°. (c) Histogram of the angular deviations: 
81% of the deviations are within the range < 9 < 30°, 96% of the deviations are within the range < 9 < 60°, and 98% of 
the deviations are within the range < 9 < 90° . 



matched titanium modified EAM (MEAM) potential^. 
A direct comparison of the errors is difficult, however, 
due to the different types of potentials and/or configura- 
tions considered in each work. 

We also determine the errors for the directions of the 
forces. The weighted average angular deviation of the 
EAM force directions from the DFT force directions is 



Axioms 

(? avg = J2 w * fl *> ( 10 ) 

i=l 

where 0i is the angle between the EAM force on atom i 
and the DFT force on atom i. Each angle in the sum is 
weighted by the corresponding scaled DFT force magni- 
tude LJi . Table Q] lists the weighted average angular devi- 
ation of the forces for each database configuration. The 
weighted average angular deviation for all the configu- 
rations is only 15.1°. The histogram in Fig. [^c) shows 
that 81% of the angular deviations are less than 30°, and 
98% of the angular deviations are less than 90°. 

The testing database contains 1,381 forces from nine 
bec configurations, one fee configuration, and one hep 
configuration. The data is generated for temperatures 
and pressures that lie between and beyond the tempera- 
tures and pressures in the fitting database. The weighted 
relative rms deviation of the EAM force magnitudes from 
the DFT values is 27%, and the weighted average angu- 
lar deviation of the EAM force directions from the DFT 
force directions is 16.5°. These values are very similar 
to the fitting database errors, indicating that the fitting 
database contains enough data to support the number of 
parameters in the potential. 



III. RESULTS AND DISCUSSION 

We assess the quality of the potential by comparing a 
wide variety of computed properties to DFT calculations 
and experimental data. All the DFT calculations use the 
same method and convergence criteria as the database 
calculations: PBE exchange-correlation functional, PAW 
pseudopotential with valence states and 4s and Ap semi- 
core states treated explicitly, T-centered fc-point meshes 
with 30 x 30 x 30 points per primitive cell, a plane-wave 
cutoff energy of 550 eV, and order-one Methfessel-Paxton 
smearing with a smearing width of 0.10 eV. We calculate 
two classes of properties: (1) properties such as elastic 
constants which are directly related to configurations in- 
cluded in the fitting database, and (2) properties such as 
surface energies which are not related to configurations 
included in the fitting database. The second class serves 
to test for over-fitting and transferability. The potential 
performs well in nearly all situations we have tested. 

A. Structural and elastic properties 

The potential's first test is reproducing the cohesive 
energy, lattice parameter, and elastic properties of bec 
Nb. We also determine the energetic stability of the bec 
lattice with respect to several other crystal structures. 
Table HIT] compares the EAM results to our DFT cal- 
culations and experimental data. The cohesive energy, 
lattice parameter, and bulk modulus are determined by 
calculating the energy of bec Nb for the volume range 
0.90Vb < V < l.lOVb, where Vq is the equilibrium vol- 
ume, and fitting the third-order Birch-Murnaghan equa- 
tion of state^r— to the results. DFT produces a co- 
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TABLE III: The EAM values for the cohesive energy, lattice 
parameter, bulk modulus, and elastic constants of bcc Nb are 
compared to DFT and experiment. The experimental lattice 
parameter and elastic constants were measured at 4.2 K. The 
EAM values for the energies and lattice parameters of the fee, 
hep, /3-W, /3-Ta, and to-Ti structures are compared to DFT 
results. The energies are relative to the energy of the bcc 
structure. 



EAM" GGA-PBE" Experiment 



Ecoh (eV/atom) 


7.09 


7.10 


7.57 6 


a (A) 


3.308 


3.309 


3.303 c 


B (GPa) 


172 


172 


173 d 


Cn (GPa) 


244 


251 


253 d 


C12 (GPa) 


136 


133 


133 d 


C 4 4 (GPa) 


32 


22 


31 d 


Ai5f cc _bcc (meV/atom) 


187 


324 




afee (A) 


4.157 


4.217 




A_Eh C p-bc C (meV/atom) 


187 


297 




dhep (A) 


2.940 


2.867 




Chcp (A) 


4.800 


5.238 




AEpw-bcc (meV/atom) 


77 


104 




CJ/3W (A) 


5.280 


5.296 




Ai^Ta-bcc (meV/atom) 


105 


83 




fl/3Ta (A) 


10.200 


10.184 




C/3Ta (A) 


5.313 


5.371 




A_E„Ti-bcc (meV/atom) 


167 


201 




a^Ti (A) 


4.845 


4.887 




C^Ti (A) 


2.735 


2.678 





"This work. 

'Experimental data from Kitted. 
Experimental data from Roberga^l. 

^Experimental data from Simmons and Wangi^. The bulk mod- 
ulus is obtained from Cn and C12: B = (Cn + 2Ci2)/3. 



hesive energy 6% lower than the experimental valued. 
The EAM cohesive energy is slightly different than the 
DFT value, since the DFT energies per atom of sev- 
eral structures under different thermodynamic conditions 
are used to construct the potential rather than the zero- 
temperature energy per atom. Both DFT and the EAM 
potential reproduce the lattice parameter measured at 
4.2 Kii with an error of less than 1%. The DFT and 
EAM bulk modulus values closely match the experimen- 
tal result^, each with an error of less than 1%. The 
bulk modulus for cubic crystals is related to the elastic 
constants Cn and C12 via B = (Cn + 2Ci2)/3. 

We compute the elastic constants C = (Cn — C\2)/2 
and C44 by straining the bcc crystal and calculating the 
resulting stress. The slopes of the stress versus stain 
curves yield the elastic constants. We use a volume- 
conserving orthorhombic strain to compute C, and a 
volume-conserving monoclinic strain for C44— . We apply 
a range of strains from —1% to +1% in each case. B and 
C determine Cn and C12. The errors of the EAM elas- 
tic constants compared to experiment^ are 4%, 2%, and 



3% for Cn, C12, and C44, respectively. The measured 
values are from single crystals at 4.2 K. In principle, the 
EAM value for C44 should closely match the DFT value 
since no experimental data is used to fit the potential. 
The POTFIT program fits to DFT stresses rather than 
the elastic constants, and the off-diagonal stress tensor 
component a xy determines C44. The off-diagonal stress 
tensor components in the fitting database are generally 
much smaller than the diagonal components, and it is dif- 
ficult to fit potentials that yield accurate C44 values. For 
example, the largest a xx value in the database is 22.4 
GPa while the largest a xy value is only 0.490 GPa. A 
large stress fitting weight must be used to produce poten- 
tials with C44 values close to the DFT and experimental 
values. 

The stability of the bcc crystal structure is demon- 
strated with respect to the fee and hep structures. DFT 
predicts that the energies per atom for fee and hep Nb 
are 323 meV and 296 meV larger than the bcc value, re- 
spectively. The EAM potential predicts that the energy 
per atom for both of these structures is 187 meV larger 
than the bcc value. Our EAM potential produces the 
ideal close-packed c/a value of 1.633 for the hep structure, 
whereas the DFT value is c/a = 1.827. For c/a = 1.633, 
the fee and hep first-nearest-neighbor distances are equal, 
as are the second-nearest-neighbor distances. Therefore, 
third-nearest-neighbor interactions must be included to 
differentiate between fee and hep for potentials with no 
angular dependence. Our EAM potential includes first- 
, second-, and third- nearest- neighbor interactions in the 
bcc structure, but only first- and second-nearest-neighbor 
interactions for the fee and hep structures. This leads to 
the ideal c/a ratio in the hep structure, and energetic 
degeneracy of the fee and hep structures. Increasing the 
range of 4> and p and including more fee and hep data in 
the fitting database produces values closer to the DFT 
results, but the bcc elastic constants, phonon dispersions, 
and vacancy formation energy agreed poorly with DFT 
and experiment. 

The bcc metals tungsten (W) and tantalum (Ta) have 
metastable f3 phases based on structures with eight atoms 
per unit cell and thirty atoms per unit cell, respectively. 
The /3-W structure has PmSn symmetry, and the /3- 
Ta structure has P^/mnm symmetry. Titanium (Ti) 
transforms from hep to bcc at 1,155 K and ambient pres- 
sure, and has a high-pressure w-phase based on a three- 
atom unit cell with P6/mmm symmetry. No data for 
the /3-W, /3-Ta, and w-Ti structures is included in the 
fitting database. The energies of these structures are 
higher than the bcc energy. The lattice parameters of all 
the structures are reproduced reasonably well with the 
largest error for the hep c value. The energetic order- 
ing of the structures is different in EAM and DFT but 
bcc is most stable in both cases. We also find that in 
finite-temperature MD simulations, the EAM potential 
stabilizes the bcc structure to the melting point for pres- 
sures below 125 GPa (see Sec. 1111 Cj) . 
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B. Point defects 

Vacancy motion is the predominant mechanism for 
solid-state diffusion, and the presence of vacancies influ- 
ences many material processes including dislocation mo- 
tion and creep. We use our EAM potential to calculate 
the single-vacancy formation energy E^ ac and migration 
energy £™ c , and the activation energy for vacancy diffu- 
sion Q vac = £^ ac + E™ c . The simulation supercells con- 
tain 8,191 atoms. We determine the vacancy migration 
energy with the nudged elastic band method*^— using 
seven image configurations between the initial and final 
configurations. The migration path is along the (111) 
direction. We also compute the vacancy energies using 
DFT for supercells with 249 atoms. In all our calcula- 
tions, the atoms are relaxed using the conjugate-gradient 
methodic. 

Table IIVI compares our EAM vacancy energies to 
our DFT results and other published EAMS.^, Finnis- 
Sinclair (F-S)51, and MEAM±i calculations. Most of the 
results are consistent with the experimental data^r— . 
Our EAM potential produces the largest formation en- 
ergy. The POTFIT program fits to the DFT energy per 
atom of each configuration in the database, instead of 
fitting to defect energies. The DFT energy-difference per 
atom between an ideal crystal and a crystal with a single 
vacancy is about 10 meV. This is close to the accuracy 
with which the EAM potential reproduces the energies 
in the fitting database and a large energy fitting weight 
must be used to achieve reasonable results. 

In the absence of strong irradiation, the equilibrium 
concentration of self-interstitial atoms in metals is much 
smaller than the concentration of vacancies. Accordingly, 
no data from configurations with interstitials is included 
in the fitting database. Instead, interstitial formation en- 
ergy calculations can test the transferability of the EAM 
potential. We determine the formation energies of six 
self-interstitial configurations: the (100) dumbbell, (110) 
dumbbell, (111) dumbbell, (111) crowdion, octahedral, 
and tetrahedral interstitials. Figure [3] shows the geome- 
try of these defects. The EAM simulation supercells con- 
tain 31,251 atoms which are relaxed using the conjugate- 
gradient method. Since no experimental data is available, 
we also compute the formation energies with DFT. The 
DFT supercells contain 251 atoms which are relaxed with 
the conjugate-gradient method. 

Table[V]lists our EAM results, along with our DFT val- 
ues and other published EAM±2, and MEAMli 
results. Our EAM potential yields self-interstitial forma- 
tion energies in the range 3.83-4.50 eV and our DFT cal- 
culations give formation energies in the range 3.95-4.89 
eV. DFT predicts that the (111) dumbbell has the lowest 
energy while our EAM potential predicts the (110) dumb- 
bell to be the lowest. The EAM results of Hu et al^ and 
the F-S results of Ackland and Thctford 5 , Rcbonato et 
al&, and Harder and BaconS also place the formation 
energy of the (110) dumbbell lowest. Since our EAM re- 
sults are not consistent with DFT, the potential may not 



(100) dumbbell 



<110> dumbbell <1 1 1> dumbbell 




FIG. 3: Schematic illustration of the (a) (100) dumbbell, (b) 
(110) dumbbell, (c) (111) dumbbell, (d) (111) crowdion, (e) 
octahedral, and (f) tetrahedral interstitials. 



be well suited for radiation damage studies. 



C. Phonon dispersion, thermal expansion, and 
pressure-volume relation 

The next group of properties relate to lattice vibra- 
tions and the thermodynamic behavior of the poten- 
tial. The calculations demonstrate the applicability of 
the potential over a large range of temperatures and 
pressures. First, we use the program phon 58 to com- 
pute the phonon spectra along high-symmetry directions 
in the Brillouin zone. The program employs the small- 
displacement method, in which atoms are moved a small 
distance from their equilibrium lattice sites. The dy- 
namical matrix obtained from the resulting forces on the 
atoms yields the phonon dispersions. 

Figure 2] compares the computed phonon spectra along 
the [£00], [£££], and [££0] directions in reciprocal space 
to experimental data- , our DFT calculations, and other 
published EAM results^. The DFT calculations are 
carried out for up to 512 atoms (8x8x8 supercells). The 
DFT results closely match the experimental data over 
much of the Brillouin zone, but the transverse modes in 
the [£00] direction show a plateau around £ = 0.25 which 
is not present in the experimental data. This discrepancy 
is not physical but rather is an artifact of the interpola- 
tion scheme used in generating the curves. The phonon 
frequencies are computed exactly at only a small num- 
ber of points in the Brillouin zone and PHON interpolates 
between these exact values to generate smooth curves. 
More exact points, i.e., even larger supercells, are re- 
quired to remove this discrepancy. Our EAM potential 
accurately describes the experimental phonon frequen- 
cies for small wave-vectors, but is unable to reproduce 
some of the features in the spectrum. This results in 
poor agreement at the zone boundaries H and N but our 
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TABLE IV: Single-vacancy formation, migration, and diffusion activation energies. The activation energy is the sum of the 
formation and migration energies: Q V ac = £™ + -E"™ c - Energies are reported in eV. 





EAM" 


GGA-PBE" 


Experiment 


EAM 6 


EAM C 


F-S d 


MEAM e 


MEAM ; 


-^vac 


3.10 


2.72 


2.6-3.1 9 


2.88 


2.76 


2.48 


2.75 


2.75 


pm 
-'-'vac 


0.77 


0.55 


0.6-1.6 9 


0.97 


0.64 


0.91 


0.54 


0.57 




3.87 


3.27 


3.6-4. l h 


3.85 


3.40 


3.39 


3.29 


3.32 



"This work. 

'EAM results of Guellil and Adams&. 

C EAM results of Hu et aiM. 

d F-S results of Harder and Bacon 51 . 

C MEAM results of Zhang et alM. 

flVlEAM results of Lee et alH. 

Experimental data from Landolt-BornsteinS. 

''Experimental data fromSi^. 



TABLE V: Formation energies for the (100) dumbbell, (110) dumbbell, (111) dumbbell, (111) crowdion, octahedral, and 
tetrahedral interstitials in eV. The lowest formation energy for each calculation is underlined. No DFT interstitial data is used 
in constructing our EAM potential. 





EAM a 


GGA-PBE a 


EAM 6 


F-S c 


F-S d 


F-S E 


MEAM / 


Moo 

-El 10 

E cr d 

Eoct 

Eld 


4.50 
3.83 
4.09 
4.02 
4.36 
4.37 


4.76 
4.31 
3.95 
3.99 
4.89 
4.56 


4.44 
4.39 
4.74 
4.93 
4.43 
4.73 


4.13 
3.99 

4.10 
4.23 
4.26 


4.821 
4.485 
4.795 
4.857 


4.85 
4.54 
4.88 
4.95 
4.91 
4.95 


2.56 



"This work. 

'EAM results of Hu et alM. 
C F-S results of Rebonato et <d&. 
d F-S results of Ackland and Thetford£. 
e F-S results of Harder and BaconS. 
^MEAM result of Lee et aZii. 



EAM results match experiment more closely than the 
EAM results of Guellil and AdamsS and Hu et alM. 

Figure Eta) shows the thermal expansion of the EAM 
potential from K to the experimental melting temper- 
ature, T^ lt = 2, 742 K. Constant-NPT MD simulations 
of 8,192 atoms at P = 1 atm yield the thermal expansion 
curve. We determine the equilibrium lattice constant 
for 138 temperatures in the range < T < 2,742 K. 
Each MD simulation runs for 500,000 steps with a 1 fs 
time step, and we determine the lattice constant for each 
temperature by averaging over the last 5,000 simulation 
steps. We compare the results to experimental data 60 
and the EAM results of Guellil and Adams£. Our EAM 
result lies just above the experimental curve while the 
Guellil and Adams potential underestimates the expan- 
sion. Our fitting database contains data for bec Nb only 
at (i) 300 K and -13 GPa to 23 GPa, (ii) 1,200 K and 2 
GPa, and (iii) 2,200 K and 7 GPa, so our results indicate 
the potential accurately interpolates to temperatures and 
pressures not included in the fit. 

Figure [S^b) shows the pressure variation in the EAM 
potential versus the relative volume V/Vq, where Vq is the 



zero-pressure volume. Constant- NPT MD simulations 
of 8,192 atoms at T = 293 K yield the pressure- volume 
curve. We determine the equilibrium volume for 50 pres- 
sures in the range < P < 125 GPa. Each MD simula- 
tion runs for 500,000 steps with a 1 fs time step, and we 
determine the volume for each pressure by averaging over 
the last 5,000 simulation steps. Zero-temperature EAM 
results are nearly identical to the 293 K values. We com- 
pare the results to data from shock experiments^ and 
our zero-temperature DFT calculations. For pressures 
to 75 GPa the agreement with DFT and experiment is 
excellent. The largest pressure in the fitting database is 
only 23 GPa from Configuration 1 in Table HI and the 
compression of bec Nb is accurately reproduced for more 
than 50 GPa beyond this pressure. The EAM result de- 
viates at larger pressures and at 125 GPa the bec crys- 
tal structure transforms to a close-packed lattice. Ex- 
periment and DFT do not show a phase transformation. 
Therefore the potential may not be well-suited for shock 
simulations, but it performs very well below 75 GPa. 
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0.0 0.5 1.0 0.5 0.0 0.5 



r Roo] h tua r [Uo] n 

FIG. 4: Phonon dispersion curves along high-symmetry directions. Our EAM results agree well with experiment^ for small 
wave vectors, but the classical EAM potential is unable to fully capture the spectrum. However, the potential improves upon 
the EAM results of Guellil and Adams 8 and Hu et al. 12 . The DFT phonon results closely match the experimental data over 
much of the Brillouin zone. The transverse modes in the [£00] direction show a plateau around £ = 0.25 which is an artifact of 
the interpolation scheme used to generate the curves. The DFT phonon results were provided by Hennig. 



D. Surface properties 

Surface properties test the transferability of our po- 
tential to configurations with low coordination-number 
since no surface data is used in constructing the poten- 
tial. Table I VII lists the relaxed surface energies for the 
{110}, {100}, and {111} surfaces. The EAM calculations 
use slab-geometry supercells with two free surfaces and 
periodic boundary conditions in the directions perpen- 
dicular to the surface normals. The conjugate-gradient 
method relaxes 600-layer slabs in directions parallel to 
the surface normals while the perpendicular dimensions 
are fixed. We compare the surface energies to our DFT 
results and published EAM 8 >^2, F-S^ 2 -, long-range empiri- 
cal potential (LREP)12, MEAM&li, and modified analyt- 
ical EAM (MAEAM)^i results. Experimental values for 
energies of individual surfaces are often based on simple 
models, so we evaluate the accuracy of the EAM surface 
energies by comparing the results to our DFT calcula- 
tions. 

We perform DFT calculation for 24-, 36-, and 48-layer 
slabs for the {100} surface, and for 12-, 18-, and 24-layer 
slabs for the {110} and {111} surfaces. A vacuum re- 
gion 10 A thick separates the periodic surface images. 
We relax the slabs in a manner similar to the EAM cal- 
culations. We use different numbers of layers to study 
the convergence of the surface energies and relaxations 
with cell size. The energy values for the different num- 
bers of layers vary by 1 meV/A 2 or less. When we in- 
crease the vacuum layer thickness to 15 A for the largest 
supercells, the surface energies change by less than 0.2 
me V/A 2 . The errors between our EAM and DFT re- 



TABLE VI: Low-index surface energies of bcc Nb in meV/A 2 
(J/m 2 ). Our EAM results closely match our DFT values, even 
though no surface data is used to construct the potential. 





£{H0} 
surf 


£{100} 
surf 


E {1 \ 1} 

SUIT 


EAM" 


127 (2.04) 


147 (2.36) 


154 (2.47) 


GGA-PBE a 


131 (2.10) 


146 (2.34) 


149 (2.39) 


EAM 6 


113 (1.81) 


123 (1.97) 




EAM C 


108 (1.73) 


120 (1.93) 




F-S d 


104 (1.67) 


122 (1.96) 




LREP e 


112 (1.79) 


131 (2.10) 


146 (2.34) 


MEAM 7 


117 (1.87) 


174 (2.79) 


126 (2.02) 


MEAM 9 


155 (2.49) 


169 (2.72) 


182 (2.92) 


MAEAM ft 


110 (1.77) 


125 (2.00) 


143 (2.28) 


"This work. 



6 EAM results of Guellil and Adams£. 

C EAM results of Hu et alH. 

d F-S results of Ackland and Finnish. 

TREP results of Dai et al—. Unrelaxed surface energies. 
■'MEAM results of BaskesS. Unrelaxed surface energies. 
fMEAM results of Lee et alH. 
ft MAEAM results of Wen and Zhangii. 



suits for the {110}, {100}, and {111} surfaces are 3.1%, 
0.7%, and 3.4%, respectively. Both methods predict that 
E^ { 0} < E^ { 0} < E^tf } . The excellent agreement be- 
tween our EAM and DFT results is surprising, consider- 
ing the fitting database does not contain configurations 
with surfaces. All the potentials give reasonable surface 
energies with respect to the DFT results. The relative 
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FIG. 5: (a) Thermal expansion curve. The thermal expansion 
of our EAM potential agrees closely with experiments^ from 
K to = 2, 742 K. Our EAM potential slightly overes- 

timates the thermal expansion, while the EAM potential of 
Guellil and Adams^ underestimates it. (b) Pressure versus 
volume curve. The experimental data are from shock experi- 
ments^. For pressures to 75 GPa, our EAM potential agrees 
well with experimental data and our DFT calculations. 



rms deviation of our EAM values from the DFT val- 
ues is 2.6%. Dai et al.'s LREP potential— has the next 
lowest relative rms deviation of 10.3%, and Lee et al.'s 
MEAM potential^ has the highest relative rms devia- 
tion of 18.9%. Baskes— unrelaxed MEAM results show a 
different ordering of the energies than DFT. 

Table IVIII lists the percent change in spacing between 
the first and second surface layers relative to the spacing 
in the bulk. We compare our EAM results to our DFT 
calculations, experimental data^, and published EAM 8 , 
F-SSa, and MEAM 11 results. Our EAM values agree very 
closely with our DFT calculations. All the methods pro- 
duce contractions of the {110}, {100}, and {111} surface 
layers, except the EAM potential of Guellil and Adams 
which predicts an expansion of the {100} layers. Our 
{100} EAM and DFT results also agree well with exper- 
iment. We do not list relaxations for layers deeper be- 
neath the surface, since the DFT results oscillate strongly 
as the number of layers in the slab changes. 



TABLE VII: Low-index surface relaxations of bcc Nb. The 
values are the relative percent change in the interplanar spac- 
ing upon relaxation. The number of layers in our slabs are in 
parentheses. 





A^ 10} (%) 


Air } (%) 


A{- } (%) 


EAM a 


-5.0(600) 


-13.9(600) 


-27.0(600) 


GGA-PF3E a 


-3.9(12) 


-12.4(24) 


-30.7(12) 


GGA-PF3E a 


-3.9(18) 


-13.0(36) 


-28.4(18) 


GGA-PBE" 


-4.5(24) 


-12.3(48) 


-27.6(24) 


Experiment 6 




-13 ±5 




EAM C 


-1.6 


+0.52 




F-S d 


-5.1 


-16.0 




MEAM 6 


-7.3 


-12.5 


-35.5 


"This work. 



''Experimental data from Lo et al£&. 
C EAM results of Guellil and AdamsIL 
d F-S results of Ackland and Finnish. 
e MEAM results of Lee et a/il. 



E. Stacking faults and dislocations 

The nonplanar core of screw dislocations in bcc transi- 
tion metals is generally accepted to be responsible for the 
complex deformation behavior of these materials^—. 
The cores of (1/2) (111) screw dislocations in bcc met- 
als spread into several planes of the (111) zone. How- 
ever, no dissociation of dislocations into well defined par- 
tial dislocations has been observed, and no metastablc 
stacking faults that could participate in such dissociation 
have been identified. The most widely used theoretical 
approach in searching for possible stacking faults is 7- 
surface calculations. The 7 surfaces represent energies 
of generalized stacking faults, formed by displacing two 
halves of a crystal relative to each other along a low- 
index crystallographic planed, i.e., the fault plane. As 
the top half of the crystal moves in the fault plane rela- 
tive to the bottom half, the crystal's ideal stacking order 
is disrupted. The resulting energy increase forms the 7 
surface, which is periodic in displacements perpendicular 
to the fault-plane normal. Minima on 7 surfaces deter- 
mine possible metastable stacking faults. 

We use our EAM potential and DFT to compute sec- 
tions through the {112} and {110} 7 surfaces in the (111) 
direction. EAM calculations with supercells containing 
60,000 atoms determine unrelaxed and relaxed 7-surface 
energies. The supercell for the {112} 7-surface has 3,000 
layers and the supercell for the {110} 7-surface has 2,000 
layers. In each case, the fault-plane divides the crystal in 
half. We calculate the energy as the top half of the crys- 
tal is displaced relative to the bottom half along (111). 
In the relaxed EAM calculations, the atoms are allowed 
to move only in the direction perpendicular to the fault- 
plane since the stacking faults are unstable. The DFT 
calculations use supercells with 24 layers for the {112} 
fault-plane and 12 layers for the {110} fault-plane to de- 
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FIG. 6: 7-surface sections in the (111) direction. The absence 
of minima indicate that there are no stable stacking faults in 
the {112} and {110} planes along this direction. The EAM 
and DFT results agree well, even though no data from configu- 
rations with stacking faults is used to construct the potential. 



termine unrelaxed 7-surface energies. 

Figure |6] compares the EAM and DFT 7-surface sec- 
tions in the (111) direction for the {112} and {110} fault 
planes. There are no minima that would indicate the 
existence of metastable stacking faults, which is con- 
sistent with 7-surface calculations for many bcc met- 
alsZ2r— . The overall agreement between our EAM and 
DFT results is very good, but the relaxed (and unre- 
laxed) EAM results show shoulders near 6/6 and 56/6 
which are absent in the DFT curves, where 6 is the mag- 
nitude of the (1/2)(111) screw dislocation Burgers vector 
b = (a/2)(lll). Relaxed 7 surfaces from F-S calcula- 
tions ' 5 show similar shoulders for the group VIB element 
Mo, but not for the group VB element Ta (Nb is also a 
group VB element). 

Figure [7J shows the two types of (1/2)(111) screw dislo- 
cation core structures found in calculations for bcc met- 
als. Figure [7J (a) shows the degenerate core, so named 
because the configurations on the left and right have the 
same energy. Figure [7J (b) shows the nondegenerate (or 
symmetric) core. Both types of cores spread into three 
{110} planes of the [111] zone. The results are presented 
using differential-displacement maps introduced by Vitek 
et alZ&. The atoms are projected onto the (111) plane, 




FIG. 7: Differential-displacement maps of the (a) degenerate 
and (b) nondegenerate (or symmetric) core-structures found 
for (1/2) [111] screw dislocations in bcc metals. For the degen- 
erate core, the structures on the left and right have the same 
energy, (c) In all cases, the core spreads into three {110} 
planes of the [111] zone. 



and the arrows represent relative atomic displacements in 
the [111] direction. The lengths of the arrows are scaled 
such that an arrow connects two atoms if its length is 
6/3. The shadings of the atoms indicate that there are 
three repeating layers of atoms in the [111] direction in 
an ideal crystal (white is the bottom layer and black is 
the top layer). 

Figure EJa) shows that our EAM potential for Nb 
produces the degenerate core. We determine the core- 
structure for a supercell containing about 900,000 atoms. 
The atoms are arranged in a cylindrical slab oriented 
such that the x axis is along the [121] direction, the y 
axis is along [101], and the z axis is along [111]. The 
radius of the cylinder is 60 nm. The supercell has 15 
(111) planes in the z direction. Periodic boundary con- 
ditions are applied in the z direction to simulate an 
infinitely-long straight screw dislocation. We insert a 
(1/2) [111] screw dislocation into the ideal crystal by dis- 
placing all the atoms in the supercell according to the 
dislocation's anisotropic elastic strain fields. The result- 
ing structure provides an initial configuration for subse- 
quent conjugate-gradient relaxation. Atoms that are less 
than 58 nm from the center of the cylinder are free to 
relax (the atomistic region) while the rest of the atoms 
are fixed at their initial positions. This fixed boundary 
condition effectively extends the dimensions of the sys- 
tem to infinity in the x and y directions. There are no 
published DFT results for the (1/2)(111) core-structure 
in Nb. F-S potentials^ produce degenerate cores for the 
group VIB metals (Cr, Mo, W) and nondegenerate cores 
for the group VB metals (V, Nb, Ta). 

Duesbery and Vitek^S. propose a criterion that relates 
the {110} 7 surface to the (1/2) (111) screw dislocation 
core structure. The criterion is based on results from F-S 
calculations and states that the degenerate core-structure 
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FIG. 8: Differential-displacement maps for the (1/2) [111] 
screw dislocation core structure produced by our EAM po- 
tential, (a) The degenerate core structure results when the 
atoms are relaxed, (b) Shear stress applied along [111] pro- 
duces a net displacement of the dislocation in the (112) plane. 
The figure shows the core-structure when the shear stress is 
just above the critical-resolved shear stress for dislocation mo- 
tion. 



forms if 



7{no} (b/3) < 27 {110} (6/6) . 



where 7{no}(V3) and 7{no}(V6) are the {110} fault 
energies at 6/3 and 6/6 along (111), respectively. Our 
EAM potential produces 7 { n }(6/3) = 0.033 eV/A 2 and 
7{no}(&/6) = 0.014 eV/A 2 . These fault energies do not 
satisfy the criterion for the degenerate core, yet this is 
the core that our EAM potential favors. This suggests 
that the Duesbery-Vitek criterion is not generally valid, 
and that the shapes of the 7 surfaces and the type of core 
structure depend on the details of atomic interactions. 

Figure |8] shows the relaxed core structure of the 
(1/2) [111] screw dislocation, and its movement under 
pure shear stress acting parallel to the Burgers vector. 




FIG. 9: Orientation of the maximum resolved shear stress 
plane (MRSSP) with respect to the (101) plane, (a) The 
angle between the MRSSP and the (101) plane is \- ( D ) Due 
to symmetry, we only need to consider —30° < x < +30° (the 
shaded interval). This range of angles includes the the (211), 
(101), and (112) planes. 



We increase the strain on the crystal in small incre- 
ments and allow the atoms in the atomistic region to 
relax after each increase in strain. The resulting shear 
stress acts in the maximum-resolved shear stress plane 
(MRSSP), and the dislocation moves when the stress 
reaches the critical-resolved shear stress (CRSS), i.e, the 
Peierls stress. We compute the CRSS for different ori- 
entations of the MRSSP. Figure ® shows that the ori- 
entations of the MRSSP are defined by the angle \ the 
MRSSP makes with the (101) plane. It is sufficient to 
consider —30° < x < +30° due to crystal symmetry. 
Figure [5Jb) shows that when the shear stress reaches the 
CRSS, the dislocation moves along the (112) plane for 
all MRSSP orientations with \ < 25°. The net mo- 
tion of the dislocation is in the (112) plane. An alter- 
native way to view this motion is the dislocation moves 
along the (101) and (011) planes in steps of (1/3) [121] 
and (1/3) [211], producing an effective slip in the (112) 
plane. The same motion is observed in atomistic sim- 
ulations of (1/2) [111] screw dislocations in Ta using a 
F-S^S potential and a model generalized pseudopotential 
theory potential^. Slip on {112} and {110} planes has 
been experimentally observed in Nb single crystals^—. 

Figure [TU] shows the CRSS for various orientations of 
the MRSSP. The results clearly demonstrate the depen- 
dence of the CRSS on the sense of shearing and illustrates 
the well-known breakdown of the Schmid law in bec met- 
als^— . This law assumes that components of the stress 
tensor other than shear in the slip plane in the slip di- 
rection play no role in the deformation process, and that 
the critical stress is independent of the sense of shearing. 
When (112) is the slip plane, the Schmid-law dependence 
of the CRSS on x has the form l/cos(% + 30°), drawn as 
a dashed curve in Fig. 1101 Deviations from the Schmid 
law becomes discernible for x ^ 15°, and rapidly increase 
as the MRSSP approaches the (211) plane. 
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FIG. 10: The critical resolved shear stress (CRSS) for disloca- 
tion motion as a function of MRSSP orientation. The figure 
shows a deviation from the Schmid law when the angle be- 
tween the MRSSP and the (101) plane is greater than about 
15°. 



F. Melting 
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FIG. 11: Two-phase melting simulations determine the melt- 
ing temperature of our EAM potential. Initially, half the sim- 
ulation cell contains liquid Nb and the other half contains bcc 
Nb. The liquid region of the simulation cell solidifies below 
the melting temperature and the solid region melts above the 
melting temperature. The figure shows the equilibrium vol- 
ume for each simulation temperature. There is a sharp jump 
in volume upon melting. 



Our primary interest is solid-state simulations, but we 
also examine melting behavior. Morris el a/.— state "for 
EAM potentials, it has been commonly observed that 
the melting temperatures are significantly lower (30% or 
more) than experimental values." Accordingly, an accu- 
rate melting temperature provides a challenging test for 
the potential. Two-phase melting simulations, in which 
the simulation cells contain solid and liquid in contact 
with each other, produce reliable melting temperatures. 
The liquid-solid interface provides nucleation sites for 
melting, thereby removing over-heating issues associated 
with single-phase melting simulations. Several methods 
based on this idea have been applied to the melting of 
metallic systems^—. 

We follow the approach of Bclonoshko et a/. 8? , in which 
constant-A^PT MD simulations determine the melting 
temperature. Initially, half the simulation cell is liquid 
and the other half is bcc. For a given pressure, we com- 
pute the average volume for a series of simulations with 
increasing temperature. The liquid region of the simu- 
lation cell solidifies below the melting temperature, and 
the solid region liquefies above the melting temperature. 
The volume of the system increases sharply across the 
melting temperature, indicating that a phase transition 
occurs. The average volume of each phase is constant 
at the melting point where the two phases coexist. Our 
simulation cells contain at least 16,500 atoms. Each MD 
simulation runs for 5,000,000 steps with a 1 fs time step, 
and we average the volumes from the last 5,000 steps. 
We check the coexistence of the phases at the melting 
temperature using at least five independent simulations. 
We find that 130,000-atom simulations produce the same 
melting temperatures as 16,500-atom simulations. 

We compute melting temperatures for simulation cells 



containing liquid in contact with a {100}, {110}, or {111} 
bcc surface. Figure QT] shows the increase in volume with 
temperature at P = 1 atm for the liquid-} 100} inter- 
face. The melting temperatures at P = 1 atm from the 
liquid-{100}, liquid-{110}, and liquid-{lll} simulations 
are 2,686 K, 2,680 K, and 2,688 K, respectively. Each 
melting temperature has an error of ±5 K. The average 
of the three melting temperatures is 2,685 K. The error 
between this value and the experimental melting temper- 
ature of 2,742 K is only 2%. The agreement is excellent 
considering that the fitting database does not contain 
data from configurations near the melting point. 

We also determine the melting curve of Nb for pres- 
sures to 2.5 GPa. Figure [T"2l shows the increase in melting 
temperature with pressure. The points are results from 
constant- NPT MD simulations, and the solid line is a 
quadratic fit through the values: T = T + aP + f3P 2 , 
where T = 2, 685.8 ± 0.2 K, a = 53.9 ± 0.3 K/GPa, and 
/3 = — 3.4±0.1 K/GPa 2 . Each data point has an error of 
±5 K. The melting curve of Nb has not been measured. 

Figure [13] shows the radial distribution functions 
(RDF) for bcc Nb at 273 K and 1 atm, and liquid Nb at 
2,750 K and 1 atm. We determine the RDFs by averag- 
ing position data from over 1,000 MD simulation steps. 
No experimental data is available for liquid Nb, so we 
compare our prediction of the liquid RDF to the result 
from an EAM potential intended for simulating liquid 
Nb^. Both potentials predict that groups of bcc peaks 
merge to form wider peaks in the liquid but there are 
small differences in the maxima of the peaks. 
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FIG. 12: The melting curve of Nb. The points are EAM 
results and the solid line fits the values. The melting tem- 
perature is computed for pressures to 2.5 GPa, with an error 
of ±5 K for each temperature. The EAM potential produces 
a melting temperature of 2,685 K at P = 1 atm. The ex- 
perimental melting temperature at this pressure is 2,742 K. 
The error between EAM and experiment is only 2%. Melting 
temperatures have not been measured for higher pressures. 



IV. 



SUMMARY 



We construct an accurate and reliable EAM potential 
for Nb as the first step in alloy potential development. 
The force-matching program potfit optimizes the EAM 
functions to a database of well-converged DFT forces, en- 
ergies, and stresses. The potential accurately reproduces 
properties tied to the fitting data, and shows excellent 
agreement with DFT and experiment for a large number 
of other quantities that are related to configurations not 
included in the fitting database. The potential describes 
structural and elastic properties, defects, and thermo- 
dynamic behavior. While the potential may not be well 
suited for shock-wave or radiation damage studies, it per- 
forms very well in all other situations we have tested. 
The potential also serves as a viable starting point for 
constructing accurate EAM potentials for Nb alloys. 
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FIG. 13: Radial distribution functions (RDF), (a) Nb is bcc 
at 293 K. (b) Liquid Nb at 2,750 K. No experimental data is 
available, so we compare the liquid result to the RDF from 
an EAM potential intended for simulating liquid Nb— . The 
two potentials produce similar results. The first and second 
neighbor shells in the solid merge into a single peak in the 
liquid, with similar behavior for higher-order neighbor shells. 



Appendix: Function modifications 

In this appendix we discuss modifications to 4>{r) and 
p(r) for small r, and to F(n) for small and large n. In 
MD simulations, fluctuations can move atoms closer to- 
gether than the minimum interatomic distance in the fit- 
ting database. The potfit program accounts for this by 
extending (j) and p to r values smaller than the inner cut- 
off radius. The cubic polynomials in the range 2.073 < 
r < 2.240 A are extended down to 1.738 < r < 2.240 A. 
Likewise, potfit extends F to n values smaller than the 
inner cutoff, and n values larger than the outer cutoff. 
The cubic polynomial in the range 0.0775 < n < 0.209 is 
extended down to —0.0186 < n < 0.209, and a very steep 
cubic polynomial is added for 1.000 < n < 1.264. Requir- 
ing continuity of F and its first and second derivatives 
at n — 1.000 determines three coefficients of the steep 
cubic function. Setting F equal to 4.828 eV at n = 1.264 
determines the final coefficient. 
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Despite these modifications, the potential is not re- 
pulsive enough for high-temperature and high-pressure 
simulations, where atoms closely approach one another. 
To overcome this limitation, we modify <f> for 1.738 < 
r < 2.073 A. We replace the function in this range by a 
steeper cubic polynomial. The continuity of <f> and its first 
and second derivatives at r = 2.073 A determines three 
coefficients. Setting the first derivative at r = 1.738 A 
equal to four times the first derivative at r = 2.073 A 
determines the final coefficient. This modified potential 
is repulsive enough at small atomic separations to pre- 
vent collapse problems for all temperature and pressure 
ranges we investigated. 



We also modify the extension of F for small n values 
to properly describe the cohesive energy. The minimum 
of the EAM energy per atom versus volume curve for bcc 
Nb equals the cohesive energy, but the embedding energy 
is not zero for n = when the atoms are far apart. There- 
fore, we replace the POTFIT modification for small n by a 
different cubic polynomial. We choose three coefficients 
to ensure continuity of the embedding function and its 
first and second derivatives at n = 0.0775. We determine 
the final coefficient by setting F(0) = 0. Table HI1 lists the 
optimized spline knots, boundary conditions, and modi- 
fications of </>, p, and F. 
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